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MODEL  VALIDATIONS  AND  PREDICTIONS  FOR  WATER  BARRIER  DEFENSE 


INTRODUCTION 

The  purpose  of  this  report  is  to  provide  both  validations  and  predictions  for  explosion  plume 
behavior.  The  experiments  presented  in  this  report  were  conducted  in  July,  1995,  in  a  water-filled 
quarry  facility  in  Rustburg,  Virginia,  operated  by  Dynamic  Testing,  Incorporated,  a  subsidiary  of  NKF 
Engineering,  Incorporated.  The  data  obtained  from  these  tests  are  used  to  validate  a  computational 
hydrodynamics  model  for  plume  predictions. 

Background 

This  report  was  prepared  in  support  of  the  Water  Barrier  Ship  Self  Defense  Concept,  managed  by 
C.  E.  Higdon  of  the  Naval  Surface  Warfare  Center,  Dahlgren  Division  (NSWCDD),  Dahlgren, 
Virginia,  Code  G23.  Under  the  sponsorship  of  the  Office  of  Naval  Research  (ONR),  Arlington, 
Virginia,  the  Center  is  developing  technology  that  has  the  potential  to  be  very  effective  in  defending 
Navy  platforms  against  high-speed,  low-flying  antiship  missiles  (ASMs).  The  concept  uses  a  “wall  of 
water”  to  provide  a  low-cost,  universal  terminal  defense  system  for  ships.  The  “wall  of  water”  or 
“water  barrier”  is  formed  from  the  shallow  detonations  of  multiple  underwater  explosives  to  protect 
the  ship  from  attacking  ASMs.  This  concept  can  be  employed  to  slow  or  stop  debris  and  warhead 
fragments  from  missiles  killed  at  very  short  range  to  preclude  significant  damage  to  the  defending 
ship.  Furthermore,  the  barrier  would  defeat  the  fusing  and  structure  of  ASMs  that  have  penetrated  the 
inner  self-defense  layer.  Further  details  of  this  concept  may  be  foimd  in  [1]. 

In  a  recent  work  [2],  observations  of  a  large  set  of  experiments  were  used  to  improve  a  set  of 
empirical  relations  for  modeling  imderwater  explosion  bubbles  in  an  incompressible  medium  (water). 
While  such  relations  were  derived  and  reported  as  long  ago  as  1948  by  Cole  [3]  these  approximations 
were  not  valid  for  very  shallow  depths  due  to  simplifying  assumptions  that  were  made.  Using  the 
Rayleigh-Plesset  equation  for  modeling  a  spherical  adiabatic  gas  bubble  oscillating  in  an  infinite 
incompressible  medium,  without  any  simplifying  assumptions,  together  with  new  measurements  of 
shallow  depth  explosion  bubbles,  the  new  empirical  relations  provide  more  accurate  initial  conditions 
for  hydrodynamics  computer  codes. 

Explosion  Dynamics 

Upon  detonation  of  an  imderwater  explosive  a  shock  wave  moving  away  from  the  charge  is 
emitted.  This  wave  reflects  off  the  surface  as  a  rarefaction  wave  which  travels  back  down  through  the 
gas  globe  of  detonation  products.  Due  to  the  tension  created  behind  the.  rarefaction  wave,  a  whitened 
area  of  cavitation  is  formed  that  rises  from  the  surface.  Under  the  surface  a  bubble  is  formed  from  the 
combustion  products,  which  expands  rapidly  due  to  the  initially  high  pressures  of  its  internal  gases. 
The  early  expansion  of  the  bubble  is  nearly  spherical,  after  which  a  water  plume  forms  above  the 
bubble.  Eventually,  the  bubble  expands  to  its  maximum  volume.  If  this  maximum  volume  has  an 
equivalent  spherical  radius  that  is  between  approximately  one  and  two  times  the  initial  charge  depth, 
a  second  jet  moving  downward  through  the  bubble  will  form  during  its  collapse  to  a  minimum  volume. 
The  duration  from  the  time  of  the  detonation  to  the  first  collapse  is  referred  to  as  the  “bubble 
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period.”  Since  the  jet  strikes  the  bottom  of  the  bubble  before  the  minimum  volume  is  attained,  the 
bubble  forms  an  annular  region.  As  the  annular  bubble  re-expands,  secondary  plumes  are  ejected 
radially,  surrounding  the  central  initial  plume. 

In  the  case  of  a  line  charge  or  several  point  charges  placed  sufficiently  close  together,  a 
cylindrical  bubble  is  formed.  The  initial  plume  forms  a  wall  of  water  above  the  line  of  charges. 
Secondaiy  plumes  erupt  on  either  side  of  the  central  plume  after  the  first  bubble  collapse. 

The  shock  related  phenomena  described  above  typically  occurs  on  the  order  of  a  tnillisecond  or 
less.  For  the  examples  discussed  in  this  report,  bubble  periods  are  approximately  0.6  s  with  secondary 
plumes  erupting  shortly  afterward.  The  entire  duration,  from  the  detonation  to  the  plume  falling 
back  down  to  the  water  surface,  usually  lasts  between  4  and  6  seconds. 

Model  Approach 

The  computer  codes  used  for  the  simulations  and  modeling  presented  in  this  report  are  based  on  a 
generalized  formulation  of  hydrodynamics  [2,4-10].  This  method  is  well  suited  for  the  study  of 
shallow-depth  explosion  plumes  for  the  following  reasons: 

1.  The  “water”  or  “liquid”  region  is  modeled  as  incompressible,  thereby  allowing  for  time  steps 
proportional  to  the  inverse  of  the  water  velocity  as  opposed  to  the  much  smaller  time  steps  that 
a  compressible  formulation  would  require  based  on  the  speed  of  sound  in  water.  This  is  important 
because  plume  behavior  occurs  on  the  order  of  seconds. 

2.  The  model  allows  for  regions  of  “spray,”  which  is  typical  of  plume  behavior  in  which  a  well 
defined  interface  between  the  bubble  and  water  or  especially  the  water  and  the  air  does  not  exist. 

3.  The  computational  model  uses  a  fixed  “Eulerian”  grid  providing  for  generality  in  studying 
complex  bubble  dynamics  and  free  surface  topology  changes.  For  shallow-depth  explosions  this 
includes  the  underwater  bubble  forming  one  or  more  annular  regions  as  a  downward  moving  jet 
intersects  the  bottom  surface  of  the  bubble  as  it  collapses,  in  addition  to  the  radial  plumes  ejected 
on  the  bubble’s  second  expansion,  and  the  eventual  venting  of  the  bubble  into  the  atmosphere. 

Our  approach  has  some  similarities  to  the  volume  of  fluid  (VOF)  approach  developed  by  Hirt  and 
Nichols  [11]  and,  more  recently,  by  Kothe  and  Mjolsness  [12].  Our  method  differs  in  some 
important  ways  from  these  typical  VOF  methods.  First,  the  theory  behind  the  model  was  designed 
specifically  for  violent  surface  motions  characterized  by  collisions  of  different  portions  of  the  free 
surface.  In  particular,  when  collisions  occur,  the  VOF  variable  will  often  attain  a  value  larger  than  one 
(due  to  numerical  error  or  fluid  elements  running  into  each  other).  When  this  occurs,  VOF 
formulations  simply  truncate  the  overage.  This  process  violates  conservation  of  mass  and  can 
introduce  small  instabilities  by  increasing  the  total  energy  of  the  solution.  Using  the  generalized 
formulation,  density  is  redistributed  in  such  a  way  that  the  total  mass  is  conserved  and  the 
momentum  is  redistributed  so  that  the  energy  is  nonincreasing  (but  may  decrease  when  liquid 
collisions  occur).  This  method  solves  conservation  of  mass  and  momentum  equations,  which  are 
subject  to  density  and  pressure  (when  cavitation  is  an  issue)  constraints.  The  density  constraint, 
together  with  the  conservation  of  mass  equation,  are  equivalent  to  the  usual  divergence  free 
constraint  for  incompressible  flow  in  regions  where  the  density  is  at  its  maximum  (liquid)  value. 
These  equations  are  solved  numerically  using  a  split  step  procedure.  First,  the  conservation  equations 
are  approximated  without  regard  to  the  constraints  using  a  second  order  Godimov  Method  with 
monotonized  slope  limiting,  as  described  in  [13].  Next,  the  density  constraint  is  imposed  through  the 
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solution  of  a  variational  inequality,  which  becomes  a  linear  complimentarity  problem  upon 
discretization.  Finally,  the  pressure  is  determined  using  a  projection  method  discretized  using  a  finite 
element  method.  This  algorithm  has  been  implemented  for  both  a  two-dimensional  (2-D)  (or  axially 
symmetric)  code  BUB2D  [2,4-7,9],  and  a  three-dimensional  (3-D)  code  in  generalized  coordinates 
BUB3D  [8-10]. 

THE  COMPUTATIONAL  MODEL 
Model  Equations 

The  computational  model  used  for  both  the  2-  and  3-D  codes  is  based  on  a  generalized 
formulation  of  hydrodynamics.  This  formulation  uses  a  fixed  spatial  domain  O,  where  the  density  p, 
velocity  u,  and  the  pressure  p  are  governed  by  the  mass  and  momentum  conservation  equations 

p,+V*(pu)  =  0  (1) 

(pu),  +  V  •  (puu)  =  -pgk  -  Vp  (2) 

subject  to  the  constraint 

P^Po.  (3) 

where  p„  is  the  constant  density  of  the  incompressible  liquid.  In  (2)  -k  is  the  unit  vector  in  the 
direction  of  the  gravitational  force,  and  g  is  the  gravitational  constant.  In  regions  where  p  =  p^, 
Eq.  (1)  becomes  the  usual  divergence  fi’ee  condition  for  incompressible  flow.  We  define  the  time 
varying  “liquid”  domain  D(/)  by 


D(0  =  {x:p(x,0  =  Po}. 


(4) 


The  non-liquid  domain  is  defined  using 

Q  -  D(0  =  A(0  u  B(0  u  C(0, 


(5) 


where  the  regions  A,  B,  and  C  are  disjoint.  Within 
uniform;  that  is 

Pa 

p{x,t)  =  \ps{t) 

Pc 


these  regions,  the  pressure  is  assumed  to  be 
X  G  A(t) 

XGB(t).  (6) 

X  G  C{t) 


In  the  above,  p^  represents  the  constant  ambient  “air”  pressure.  The  “bubble”  pressure,  Pgit),  is 
usually  determined  using  an  adiabatic  gas  assumption 

PB(t)=c(V,(t)y;  (7) 

where  c  is  constant,  7  is  the  (constant)  ratio  of  specific  heats  of  the  bubble  gases,  and  Vg(t)  is  the 
bubble  volume,  which  can  be  determined  using 


Finally,  Pc  is  the  “cavitation”  pressure,  which  is  usually  set  to  the  vapor  pressure  of  the  liquid  at 
some  specified  temperature.  When  cavitation  is  to  be  modeled,  an  additional  constraint  is  imposed  on 
the  pressure,  namely  that  p(x,t)>  Pc  (see  [8]).  For  the  results  presented  here,  this  constraint  was 
not  imposed,  that  is,  Pc  =  -°°- 

Numerical  Algorithm 

Assume  that  the  density  and  velocity,  p",u"  at  time  step  n  are  known  together  with  the  pr^sme 
gradient  at  the  previous  half  step,  .  This  solution  is  evolved  from  time  t  =  t  —>  t  +T  =  t 

using  the  following  three  step  time  split  procedure. 


Convection 

The  solution  is  first  advanced  (p",u")->  (p,u)  by  “solving”  the  conservation  laws  (1-2)  without 
including  the  term  Vp  on  the  right-hand  side  of  Eq.  (2)  and  without  regard  to  the  constraint  Eq.  (3). 
This  step  is  fully  discretized  using  a  formally  second-order  Godunov-type  method,  which  uses  slope 
limiting  in  space  and  explicit  predictor-corrector  time  stepping  (e.g.,  [13]).  Although  the  pressure 
gradient  is  not  explicitly  added  to  the  momentum  here,  it  is  used  within  the  predictor  step  of  the 
Godimov  method.  Further  details  of  this  step  may  be  found  in  Refs.  4  and  6  for  axially  symmetric 
flow  problems. 

Redistribution  of  Density  and  Momenta 

Next,  the  density  and  momenta  are  redistributed  ^p, uj— >  (p, u  )  so  that  the  constraint  Eq.  (3)  is 

satisfied,  the  global  conservation  of  mass  and  momenta  are  maintained,  and  the  energy  is 
nonincreasing.  The  density  is  redistributed  using  an  approximate  solution  to  the  obstacle  problem 


jPo-P 

1  0 


if  .^>0 
if  i/  =  0’ 


(9a) 


and  setting 

p  =  p  +  V^H. 


(9b) 


These  equations  have  been  derived  by  considering  a  Boltzmann  formulation  for  modeling  inelastic 
fluid  collisions  and  are  discussed  in  greater  detail  in  Refs.  4,  6,  and  9.  These  references  also  contain 
details  of  the  numerical  discretization  of  Eq.  (9)  and  the  solution  procedure,  which  eniploys  a 
constrained  direction  preconditioned  conjugate  gradient  method.  The  momenta  redistributions  are 
determined  as  solutions  of  two  (or  three  for  3-D  problems)  elliptic  self-adjoint  problems 

^  =  ^  +  V^(//5).  (10) 


Discretizations  of  Eq.  (10)  yield  systems  of  linear  equations  with  diagonally  dominant  matrices, 
which  are  efficiently  solved  using  a  diagonally  preconditioned  conjugate  gradient  method.  The 
importance  of  this  step  to  the  overall  accuracy  and  stability  of  the  algorithm  was  discussed  in  Ref.  7. 
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As  mentioned  in  the  introduction,  this  redistribution  step  is  a  major  distinguishing  feature  between 
this  algorithm  and  other  VOF  approaches  (e.g.,  [11,12])  which  would  simply  truncate  the  density 
values. 

Multiple  Bubbles  and  Venting 

After  the  redistribution,  p  =  and  the  new  nonliquid  regions  are  then  determined  along  with 
the  pressure  in  each  of  its  connected  subsets.  In  the  computational  space,  the  new  liquid  region, 
D"^'  =  D(r'^' ),  is  defined  to  be  the  collection  of  grid  cells  C,  such  that 

^(l-£p)Po.  (11) 

where  pj^'  is  the  density  in  cell  C,.  For  shallow-depth  explosion  bubbles,  the  choice  of  is 
important.  In  general,  small  values  of  will  cause  cells  with  only  slightly  less  density  than  the  liquid 

to  be  treated  as  regions  with  uniform  pressure,  while  larger  values  will  cause  more  of  the  “spray” 
(where  0<  p  < Po)  to  be  treated  as  a  variable  density  incompressible  region. 

In  addition  to  the  distinction  in  the  nonliquid  regions  designated  by  “air”  A,  “bubble”  B,  and 
“cavitation”  C,  each  component  (connected  disjoint  subset)  of  B  is  also  accoimted  for.  That  is, 
suppose 


K" 

Bn  _  nn 

h^\ 

where  K”  is  the  number  of  distinct  bubbles  at  time  step  n,  and  is  the  component  corresponding  to 
bubble  k  at  time  step  n.  When  the  bubble  components  remain  distinct,  their  pressures  behave 
adiabatically 


pl-c,{y,y, 

where  represents  the  volume  of  and  is  constant  for  all  steps  n  for  which  the  component 
has  no  interactions.  If  a  bubble  component  splits  into  two  distinct  regions,  say,  B^  ->  5,"^'  the 

new  pressures  and  p^'  are  computed  assuming  the  volume  changes  occurred  before  the  split; 
that  is, 

Pi  Pm 

Similarly,  if  two  distinct  bubbles  merge  into  a  single  component  region,  for  example,  B"uB^  5^*, 
the  new  pressure  is  given  by 

(P"v," + p:v:wr +v:y-^ 

'  (vry 

These  formulas  have  been  extended  to  the  general  case  treating  any  finite  number  of  bubbles 
merging  and  splitting  in  [5].  Whenever  merging  occurs,  the  pressure  of  the  new  component  changes 
instantaneously.  Similarly,  when  a  bubble  comes  in  eontact  with  the  air  region  (that  is  a  cell  in  B^  is 
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adjacent  to  a  cell  in  A),  the  bubble  is  said  to  “vent”  into  the  air  with  the  pressure  instantaneously 
changing  to  the  air  pressure.  This  is  obviously  a  crude  approximation  to  the  finite  amount  of  time 
the  venting  would  actually  take  or  to  the  partial  venting  of  the  bubble.  In  particular,  if  the  relatively 
thin  layer  of  water  between  the  bubble  and  air  is  under-resolved  by  the  computational  grid,  the 
computed  bubble  may  vent  prematurely  causing  gross  errors  in  the  subsequently  computed  dynamics. 
In  order  to  prevent  this  premature  venting,  the  value  of  Sp  in  (1 1)  is  allowed  to  depend  upon  both 
time  and  space  by  the  prescription 

QeN(A")^  (12) 

^  ^  [eg  otherwise 

where  N(A")  is  the  set  of  cells  that  are  either  in  or  adjacent  to  the  time  varying  domain  A". 

Consider  the  case  where  p;>(l-er’)po  so  that  C,gD”  but,  on  the  next  step,  the  density 

decreases  below  the  liquid  cutoff,  that  is,  pr’  <(l-e,”)po-  If  there  is  a  cell  in  the  neighborhood  that 

was  nonliquid  (C*  GN(Q)n(A"  uB"  uC"))  then  cell  C,  would  be  merged  into  the  nonliquid 

component  eontaining  that  cell.  If  there  are  more  than  one  nonliquid  components  in  the 
neighborhood,  these  components  are  merged  together.  However,  the  case  when  there  are  no 
nonliquid  neighboring  cells  is  more  problematic.  Since  the  velocity  filed  in  the  liquid  is  divergence 
free,  the  density  can  decrease  only  because  of  numerical  error.  Previously  this  case  was  treated  by 
resetting  p,"^'  =[l-£"jpo  so  that  the  eell  remained  a  liquid  cell.  However,  this  obviously  adds  mass, 

and  violates  the  conservation  law.  Furthermore,  this  added  mass  can  be  large  in  the  special  case  when 
=  but  e”=£g.  Even  without  a  decrease  in  density,  the  added  mass  in  this  case  would  be 
-£g)po  IQ  I,  where  |Q  |  is  the  volume  of  Q.  In  order  to  avoid  these  potentially  large  errors  in 

added  mass,  the  density  is  no  longer  modified  in  this  case.  Instead,  the  cell  remains  part  of  the  liquid 
domain,  Q  e  .  Thus,  the  liquid  domain  is  now  defined  as  the  collection  of  cells  Q  such  that 

either  Eq.  (11)  holds  or  N(Q)cD”.  That  is: 

D'^'=UQ:pr‘>(l-£p)po  or  N(Q)cD”.  (13) 


Pressure  Projection 

In  the  nonliquid  region  u  =  .  However,  u  is  not  consistent  with  Eq.  (2)  in  the  new  liquid 

region.  In  this  region  the  velocity  is  corrected,  using  the  gradient  of  the  new  pressure, 

V^n+i/2  pressure  /?”^^^^is  the  solution  of 


yp"') 


=  V  •  u  in  D 


(14) 


The  new  velocity,  given  by 


n 


n+l 


U-T 


Vp^ 


(15) 
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is  divergence  free  in  D"^'  and,  thus,  is  consistent  with  Eq.  (1).  Equations  (14)  and  (15)  define  a 
projection  u"^’  =  P(u)  onto  the  space  of  divergence  free  velocities.  This  equation  is  discretized  using 
a  finite  element  method  with  bilinear  (in  2-D)  or  trilinear  (in  3-D)  elements.  This  spatial 
discretization  produces  an  approximate  projection  which  has  been  analyzed  in  Ref.  14.  The  resulting 
linear  system  from  the  discretization  of  Eq.  (14)  is  solved  using  an  incomplete  Cholesky 
preconditioned  conjugate  gradient  method. 

To  determine  the  pressure  uniquely  using  Eq.  (14),  boundary  conditions  must  be  specified.  On 
those  portions  of  the  boundary  of  that  correspond  to  “wall”  boimdaries  of  O,  a  Neumann 
condition  is  specified,  namely 
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Note  that  this  condition,  together  with  Eq.  (15),  implies  that  u"‘^’«n  =  0  along  wall  boundaries.  Since 
the  pressure  is  assumed  to  be  continuous,  Dirichlet  conditions  for  the  pressure  along  the  nonliquid 
regions  are  determined  according  to  Eq.  (6).  In  particular,  we  specify 
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along  those  boimdaries  common  to  the  nonliquid  regions.  Hydrostatic  pressure  Dirichlet  conditions 
are  set  on  other  portions  of  the  boundary  of  Q  to  model  “cutoff’  boimdaries. 

Initial  Conditions 

For  a  single  charge,  the  bubble  is  initialized  as  a  spherical  “void,”  with  zero  density,  radius  A^, 
and  pressure  pi,  surrounded  by  a  liquid  region  at  rest.  In  the  case  of  a  line  charge  or  a  series  of  single 
charges  placed  sufficiently  close  in  a  straight  line,  the  bubble  is  initialized  as  a  circular  cylinder  of 
infinite  length.  The  initial  values  for  the  bubble  radius  and  pressure  depend  on  the  hydrostatic 
pressure  at  the  depth  of  the  charge,  the  charge  weight,  and  empirical  constants  which  depend  on  the 
charge  type,  derived  from  considering  the  equation  of  motion  of  a  spherical  bubble  in  an  infinite 
incompressible  medium.  In  this  case,  the  bubble  remains  spherical,  and  its  radius  oscillates  periodically 
between  its  minimum,  and  maximum,  A^,  values. 

Empirical  Relations 

The  following  empirical  relations  are  valid  for  shallow-depth  underwater  explosions  (see  Ref.  2) 
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pl=ps^-Yi 
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In  the  equations  listed  above,  q  is  the  empirical  minimum  radius  constant,  which  depends  on  the 
charge  type;  W  is  the  charge  mass;  d  is  the  initial  charge  depth;  and  is  the  ambient  air  pressure. 
If  d  is  measured  in  units  of  feet  and  in  feet  of  water,  then  is  the  hydrostatic  pressure  at  the 
charge  depth  in  units  of  feet  of  water.  In  Eq.  (20),  y  is  the  ratio  of  specific  heats  for  the  bubble 
gases,  and  is  the  empirical  maximum  radius  constant.  In  the  above  equations,  the  value  for  a 
must  be  determined  as  the  solution  to  Eq.  (20)  with  a  given  value  for  the  right-hand  side.  This  can  be 
done  approximately  using  Newton’s  method.  The  value  cc  is  also  the  ratio  of  maximum  to  minimum 
bubble  radii, 

a  =  ^,  (22) 
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In  cases  when  the  depth  ranges  between  100  and  1000  ft,  it  was  noted  in  Ref.  3  that  values  for 
remain  nearly  constant.  Because  of  this  and  the  interest  in  the  relatively  deeper  charge  depths, 
subsequent  reports  (e.g.,  [15,16])  used  Eq.  (23)  under  this  assumption.  However,  depends  on  a, 
which  can  change  significantly,  particularly  at  shallow  depths.  This  dependence  was  studied  in  detail 
in  Ref.  2. 

For  line  charges  approximated  as  a  circular  cylinder  of  infinite  extent,  the  corresponding 
formulas  were  derived  in  Ref.  2. 
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In  Eqs.  (24)  and  (28),  M  is  the  mass  per  unit  length. 

For  the  tests  considered  in  this  study,  all  charges  were  Composition  C-4.  Based  on  our  previous 
analysis  (Ref.  2,  Table  3-7),  the  values  used  in  this  study  are  y  =  1.34,  qr  =  0.286 ,  and  =15.3. 
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Shock  Effects 

It  was  noted  by  Kedrinskii  in  Ref.  17  that  using  an  incompressible  liquid  model  for  shallow-depth 
explosion  simulations  generally  underpredicts  the  plume  heights  for  early  times.  It  has  been 
demonstrated  in  Refs.  2  and  17,  that  an  indentation  of  the  free  surface  directly  above  the  charge  will 
increase  the  plume  heights  predicted  by  an  incompressible  liquid  model.  This  indentation  represents 
the  effects  of  spalling  from  the  reflection  of  the  detonation  shock  wave  from  the  surface.  In  addition 
to  this  reflection,  shock  interaction  from  the  simultaneous  detonation  of  multiple  charges  has  been 
shown  to  cause  plume  fingering  between  the  initial  charge  locations.  An  empirical  model  for  this 
phenomenon  is  shown  in  Fig.  1. 


NUMERICAL  RESULTS  AND  VALIDATIONS 

In  this  chapter,  validations  of  the  2-D  computational  model  are  first  presented.  Comparisons  of 
the  computations  to  experimental  data  include  plume  height  measurements  from  video  cameras 
images  and  plume  density  measurements  from  both  conductivity  probes  and  microwave  absorption. 
The  computational  model  is  then  used  to  determine  an  “optimal  depth  study”  in  which  a  measure  of 
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the  plume  density  as  a  function  of  charge  depth  is  maximized.  Finally,  the  empirical  shock  model 
displayed  in  Fig.  1  is  validated  using  the  BUB3D  code. 

The  tests  described  here  were  conducted  in  a  130-foot-deep  water-filled  quarry  operated  by 
Dynamic  Testing,  Incorporated  (DTI)  in  Rustburg,  Virginia,  during  July  1995.  The  charges  used  in 
the  experiments  were  configured  in  a  line  comprised  of  five  to  eight  discrete  10-lb  blocks  of  C-4 
separated  by  8  ft  or  as  a  continuous  line  charge  of  length  40  ft  to  56  ft.  A  more  detailed  description 
of  these  tests  has  appeared  in  Ref.  18.  The  subset  of  tests  considered  here  are  summarized  in  Table  1. 
In  this  table,  N  is  the  number  of  discrete  charges  in  the  configurations,  S  is  the  distance  between  the 
center  of  the  charges,  is  the  maximum  theoretical  radius  of  the  cylindrical  bubble  as  determined 
using  Eq.  (28),  and  C  is  the  scaled  depth,  C  =  dl  The  parentheses  surrounding  the  number  of 
charges  of  the  continuous  type  indicate  that  these  were  actually  comprised  of  individual  1-lb  blocks 
of  C-4  in  a  line. 


Table  1  —  Test  Shot  Descriptions 


Shot 

No. 

Type 

^(ft) 

N 

^(ft) 

M 

(Ib/ft) 

A2D) 

^max 

(ft). 

C 

2 

6 

7 

9 

10 

11 

Discrete 

Discrete 

Continuous 

Discrete 

Discrete 

Continuous 

8.2 

8.2 

8.2 

8.2 

9.4 

9.4 

5 

8 

(56) 

5 

5 

1401 

8 

8 

(1) 

8 

8 

(1) 

1.25 

1.25 

1.25 

1.25 

1.25 

1.25 

11.53 

11.53 

11.53 

11.53 

11.36 

11.36 

0.71 

0.71 

0.71 

0.71 

0.83 

0.83 

Validations  of  the  Two-Dimensional  Model 

The  two-dimensional  (2-D)  model  uses  an  initial  approximation  of  the  bubble  as  a  circular 
cylinder  of  infinite  extent.  This  approximation  is  reasonable  for  the  discrete  tests  only  if  the 
distance  between  charges  is  relatively  small  when  compared  to  the  maximum  diameter  of  the 
spherical  bubble  from  an  individual  charge.  For  the  discrete  tests  listed  in  Table  1,  the  individual  10- 
pound  charges  of  C-4  would  generate  bubbles  with  maximum  radii  A„^  =  9.3,  when  d-%2,  and 
A  =9.2  at  d  =  9A.  Since  these  values  are  greater  than  the  charge  standoff  of  8  ft,  it  can  be 

ftlCDC  ^  , 

expected  that  the  bubble  dynamics  will  be  adequately  represented  by  the  2-D  approximation. 
However,  it  can  be  expected  that  shot  numbers  7  and  11  will  be  better  represented  by  this 
approximation. 

The  grids  used  in  this  study  are  tensor  product  grids.  That  is,  the  grid  point  locations  may  be 
defined  as  a  tensor  product  of  two  one-dimensional  (1-D)  grids.  Typically,  the  grids  are  constructed 
with  uniform  spacing  in  both  directions  in  the  vicinity  of  the  bubble,  with  grid  stretching  to 
approximate  either  conditions  at  infinity  or  a  wall  boundary  relatively  far  away.  The  x-direction  is 
taken  to  be  the  horizontal  line  perpendicular  to  the  line  charge  (cylindrical  axis).  For  single  line 
charges,  x  =  0  is  a  symmetry  plane.  The  z  coordinate  measures  vertical  displacement  with  the 
convention  that  z  =  0  corresponds  to  the  initial  location  of  the  water-air  interface. 

Cylindrical  explosion  bubbles  are  more  difficult  to  resolve  numerically  than  axially  symmetric 
bubbles  due  to  the  greater  values  of  the  radius  ratio  For  the  test  shots  in  Table  1,  the  initial 
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(minimum)  bubble  radius  (determined  using  Eq.  (24))  is  =  0.197  so  that  >  57.  Thus,  the  use 
of  a  single  uniform  grid  capable  of  resolving  the  initial  bubble  and  extending  beyond  the  bubble  at  its 
maximum  size  would  contain  a  prohibitively  large  number  of  cells.  One  method  of  alleviating  this 
problem  is  through  the  use  of  two  separate  grids. 

The  2-D  solutions  were  computed  in  two  steps.  Initially,  a  grid  that  was  fine  in  a  region 
surrounding  the  charge  line  was  used  imtil  the  bubble  approached  the  boimdary  of  the  fine  region. 
Then  the  solution  was  remapped,  conserving  mass  and  momentum  onto  a  grid  that  was  coarser  than 
the  fine  grid  region  of  the  first  grid  but  still  able  to  resolve  the  bubble  after  the  initial  grid  had  been 
used.  This  second  grid  was  uniform  in  a  large  enough  region  to  contain  the  important  long-time 
dynamics  of  the  problem.  For  example,  the  initial  grid  used  in  these  computations  consisted  of 
40  X  100  square  cells  of  size  =0.05  in  the  region  0<jc<2  intersected  with  -d -2.5  <  z  <-d +2.5 . 
Outside  of  the  uniform  region,  cells  were  stretched  horizontally  to  x=X^  =  ll0  and  downward  to 
z  =  Zj  =  -115  using  40  additional  grid  lines  in  each  direction.  Above  the  uniform  region,  the  spacing 
in  the  z  direction  wasxmiformly  set  to  Az  =  0.1  in  the  region  -<7+2.5  <z  <2.  The  initial  grid  was 
used  for  0  <  t  <  t,.  =  0.007,  while  the  cylindrical  bubble  grew  from  its  initial  radius  of  0.197  to 
approximately  1 .6.  The  computed  solution  at  t  =  ti  was  then  remapped  onto  the  second  grid. 

The  numerical  errors  were  approximated  by  using  two  different  grids  for  the  second  step  of  the 
computations.  Each  grid  contained  square  cells  of  size  in  the  region  0  <  x  <  12  intersected  with 
-<7-11  <z <28.  As  with  the  initial  grid  domain,  the  second  grids  were  stretched  horizontally  to 
x  =  X^  =  110  and  downward  to  z  =  Zj  =-115.  To  simulate  the  tests  where  the  initial  charge  was  at 
depth  <7  =  8.2,  the  grids  were  stretched  upward  to  z  =  Z,  =180  and,  in  the  case  when  <7  =  9.4,  they 
were  stretched  to  z  =  Z,  =125.  The  “fine”  grids  contained  a  imiform  region  with  =  0.2  and  used  a 
total  of  100  X  305  cells,  and  the  “coarse”  grids  contained  a  uniform  region  with  =  0.4  and  used  a 
total  of  50  X  160  cells. 

The  effect  of  using  an  indented  channel  in  the  surface  above  the  line  of  charges  was  also  tested. 
This  was  done  by  performing  the  computations  using  both  an  initially  “flat”  and  “indented”  free 
surface.  The  indentations  are  determined  using  the  model  shown  in  Fig.  1  with  Rj=l  .03d,  which  is 
the  same  value  used  in  Ref.  2  for  axially  symmetric  computations  with  a  single  charge. 

A  summary  of  the  runs  with  computed  bubble  periods,  7),  (the  computed  time  that  the  bubble 
attains  its  minimum  radius)  and  maximum  radii,  is  listed  in  Table  2.  Here,  the  radius  refers  to 

the  equivalent  radius  of  a  cylinder  with  the  same  cross-sectional  area  since  the  bubbles  will  in  general 
not  remain  circular.  The  number  of  time  steps  taken  to  the  first  bubble  minimum,  N(T^),  and  the 
total  number  of  steps,  N{Tp),  to  reach  the  final  time  Tj,  =6.0,  are  also  listed.  The  time  steps  for  the 
runs  were  adaptively  selected,  based  on  changes  in  equivalent  radius.  The  number  of  steps  taken  until 
the  first  bubble  minimum  is  determined  by  an  input  tolerance,  which  was  halved  for  the  fine  grid  runs. 
The  execution  times  in  Table  2  are  in  minutes  on  a  Silicon  Graphics  Impact  2  workstation  with  an 
RIOOOO  processor,  compiled  using  y77  -03  -mips4  -n32  -r 10000.  The  times  listed  are  for  the 
completion  of  the  second  step  of  the  run  only.  The  execution  time  for  the  initial  grid  was 
approximately  15  min. 
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Table  2  —  Run  Summary 


Run 

d 

Surface 

Second  Grid 

J(2D) 

T, 

MW 

MW 

Execution 

time 

1 

flat 

indented 

flat 

indented 

coarse 

coarse 

fine 

fine 

10.824 

10.761 

10.903 

10.820 

0.5921 

0.5898 

0.5862 

0.5859 

466 

448 

860 

874 

1912 

1978 

3574 

3605 

11-C 

11 -Cl 
11-F 

11 -FI 

EQ 

mKm 

Wmm 

flat 

indented 

flat 

indented 

coarse 

coarse 

fine 

fine 

0.6087 

0.6069 

0.6034 

0.6019 

463 

465 

876 

884 

1927 

1955 

3826 

3761 

28.4 
28.7 

402.4 
415.8 

The  values  of  listed  in  Table  2  indicate  that  the  depth  has  little  influence  on  the  maximum 

radius.  However,  the  fine  grid  values  are  approximately  1%  greater  than  the  coarse  grid  values.  Since 
the  grid  sizes  and  time  steps  were  approximately  halved  in  the  fine  grid  and  first  order  convergence 
has  been  observed  for  the  maximum  radius  (and  period)  using  this  numerical  method,  the  difference  in 
the  two  solutions  provides  an  estimate  of  the  numerical  error.  The  computed  periods  at  the  shallower 
depth  are  consistently  shorter  than  at  the  greater  depth.  This  is  consistent  with  the  theoretical  and 
experimental  observations  reported  in  Ref.  2  due  to  the  proximity  of  the  free  surface.  As  with  the 
maximum  radius  values,  the  periods  computed  using  the  fine  grids  are  approximately  1%  greater  than 
on  the  corresponding  coarse  grid 

Plume  Observations 

Density  contours  of  the  plumes  at  t  =  0.2,  0.5,  1.0,  and  2.5  for  the  computational  Runs  7-C,  7-F, 
7-CI,  and  7-FI  are  displayed  in  Fig.  2.  These  images  are  based  on  a  shading  corresponding  to  the 
logarithm  of  the  density.  This  enables  the  relatively  low  density  of  the  spray  in  the  plume  to  be 
easily  observed.  Densities  less  than  0.1%  of  the  water  density  are  represented  as  a  “white”  region,  and 
“black”  corresponds  to  the  water  density. 

At  time  t  =  0.2,  the  “egg-shaped”  bubble  has  nearly  attained  its  maximum  radius.  The  bubble 
profiles  for  the  four  cases  at  this  time  are  nearly  indistinguishable,  except  for  the  slightly  irregular 
profile  of  the  bubble  in  Run  7-C.  A  small  region  of  “numerical  spray”  has  begun  to  form  above  the 
water  on  top  of  the  bubble.  This  corresponds  to  the  “Rayleigh-Taylor”  instability  at  the  water-air 
surface,  where  the  water  has  an  upward  velocity  but  is  being  accelerated  downward  due  to  the  below 
ambient  pressure  inside  the  bubble.  Further  discussions  of  instabilities  due  to  explosion  bubble 
dynamics  are  included  in  Refs.  4  and  9.  A  small  thin  upward  moving  water  jet  has  formed  in  the  cases 
when  the  initial  surface  was  indented.  This  phenomenon  has  been  described  in  Ref  17  and  examined 
further  in  Ref.  2. 

At  time  t  =  0.5,  the  bubble  has  begun  to  collapse.  During  this  collapse  a  high-pressure  stagnation 
point  forms  in  the  water  above  the  bubble.  This  in  turn  causes  the  formation  of  a  “double  jet”;  a 
small  one  moving  downward  through  the  bubble,  and  a  larger  one  forming  the  central  water  column  of 
the  plume.  Similar  jetting  behavior  has  been  observed  for  axially  symmetric  bubbles  by  Blake  and 
Gibson  in  Ref  19.  In  the  computations,  this  central  column  is  surrounded  by  the  numerical  spray 
region,  which  is  moving  upward  but  is  being  accelerated  downward  by  gravity.  The  thin  jet  caused  by 
the  indentation  rises  approximately  10  ft  above  the  numerical  spray.  Once  again,  there  is  little 
difference  between  the  bubble  and  plume  profiles  for  the  coarse  and  fine  grid  run  comparisons. 
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reached  a  height  of  over  25  ft,  and  a  total  width  of  over  40  ft.  At  this  time,  differences  in  the  details 
inside  the  plume  structure  become  more  evident  at  the  different  grid  resolutions.  On  the  fine  grid 
runs,  the  plume  structures  are  almost  identical,  except  for  the  thin  central  plume,  which  remains 
higher  in  the  indented  surface  case.  The  top  of  the  secondary  plume  from  Run  7-C  is  slightly  wider 
and  lower  than  with  Run  7-CI.  Differences  in  the  solutions  as  time  progresses  are  expected  due  to  the 
water-air  instabilities  mentioned  above  and  the  instabilities  at  the  bubble  interface  near  bubble 
minimums.  However,  the  overall  qualitative  behavior  of  the  plume  and  bubble  appears  to  be 
reproducible. 

Plume  differences  are  more  apparent  at  t  =  2.5.  At  this  time,  the  plumes  have  begun  falling 
downward  with  gravity,  and  the  bubble  has  little  remaining  energy.  Bubble  energy  is  reduced  in  our 
computational  model  through  both  numerical  dissipation  and  through  our  redistribution  step  and 
treatment  of  liquid-on-liquid  collisions.  For  further  discussions  on  energy  losses,  see  Refs.  4,  6,  7, 
and  9.  The  widths  and  heights  of  the  secondary  plumes  for  the  coarse  grid  runs  have  continued  to 
diverge.  In  the  indented  runs,  the  plume  structure  is  similar  between  the  two  grids,  with  the  coarse  grid 
secondary  plume  extending  higher  (probably  due  to  the  use  of  the  stretched  grid  cells).  Some  detail  in 
the  plume  structure  at  a  height  of  about  30  ft  appears  in  the  fine  grid  Run  7-FI  due  to  one  of  the 
bubble  pulsations.  During  the  computations,  the  bubble  typically  underwent  ten  pulsations  before 
venting  at  about  t  =  4.5.  This  plume  detail  can  also  be  seen  at  a  height  of  about  20  ft  in  Run  7-F  ^d 
just  above  the  surface  for  the  coarse  grid  runs.  Comparing  the  four  runs  demonstrates  that  the  major 
difference  between  using  the  flat  or  indented  surface  is  the  height  of  the  thin  central  plume,  and  even 
though  some  detail  is  lost  with  the  coarse  grids,  their  solutions  agree  qualitatively  with  the  finer  grid 
runs,  particularly  at  the  early  times. 

Figure  3  displays  graphs  of  the  computed  plume  heights  as  a  function  of  time  for  DTI  Shot  7. 
Here,  the  plume  height  was  defined  as  the  highest  location  of  a  grid  cell  having  a  density  greater  than 
1%  of  the  water  density.  When  the  free  surface  has  an  initial  indentation  the  primary  central  plume 
always  remains  higher  than  the  secondary  plumes,  as  indicated  by  the  smooth  nearly  parabolic 
profiles  fi-om  Runs  7-CI  and  7-FI.  In  this  case  the  difference  between  the  fine  and  coarse  grids  is 
relatively  small.  For  the  first  second  the  relative  error  is  less  than  3%.  The  error  increases  to  about 
6%  at  t  =  2  s,  and  the  error  in  the  maximum  computed  plume  height  is  under  8%.  The  inflection 
point  in  the  plume  height  for  Run  7-F  shortly  after  t=  1.1  was  caused  by  the  secondary  plume  rising 
above  the  central  plume.  In  Fig.  2  the  details  of  this  secondary  plume,  and  of  subsequently  ejected 
plumes  are  not  reproduced  on  the  coarse  grid  Run  7-C.  The  discontinuity  in  the  plume  height  at 
t  =  3.5  s  for  Run  7-F  occurs  as  the  density  of  the  falling  plume  decreases  below  1%  of  the  water 
density.  This  happens  because  the  plume  has  a  horizontal  velocity  and  becomes  under-resolved 
(diffused)  as  it  passes  through  the  stretched  cells  in  the  grid. 

Figure  4  displays  the  computed  plume  heights  for  Shot  11.  The  plume  heights  are  roughly  25 /o 
lower  than  the  computed  results  for  Shot  7.  As  before  the  agreement  between  the  coarse  and  fine  grid 
results  with  the  initially  indented  free  surface  (Runs  11-CI  and  11-FI)  is  very  good,  except  for  t>  3, 
when  part  of  the  coarse  grid  secondary  plume  extends  above  the  primary  plume. 

Measurements  of  the  plume  heights  using  video  frames  were  documented  in  Ref.  18.  The 
measured  heights  for  Shot  7  are  compared  to  the  results  from  the  fine  grid  computations  in  Fig.  5. 
After  t  =  1.2  s,  plume  height  measurements  were  not  possible  due  to  a  combination  of  the  poor 
definition  of  the  top  of  the  plume  and  low  contrast.  While  the  computations  using  a  flat  initial 
surface  under  predict  the  measurements  by  over  30%,  the  use  of  the  initial  indentation  produces 
relative  errors  under  10%.  At  early  times  the  higher  values  for  the  measured  heights  are  conjectured 
to  be  caused  by  the  initial  shock  spalling  water  from  the  surface.  Recall,  we  are  not  modeling  the 


Model  Validations  and  Predictions 


15 


shock  directly,  but  only  empirically  using  the  indentation  above  the  surface.  At  later  times,  the 
computed  heights  for  Rim  11 -FI  overtake  the  measurements.  This  is  very  likely  due  to  drag  on  the 
top  of  the  plume  in  the  air;  another  phenomenon  that  is  not  included  in  our  model. 


Time  (s) 


Fig.  3  —  Computed  plume  heights  for  Shot  7 


Time  (S) 


Fig.  4  —  Computed  plume  heights  for  Shot  11 


The  results  for  the  deeper  case.  Shot  11,  are  shown  in  Fig.  6.  As  before,  the  computation  with 
the  flat  initial  surface.  Run  11-F,  under  predicts  the  measured  values  at  all  times,  while  Run  11 -FI 
predicts  heights  below  the  measurements  at  early  times  and  exceeding  the  measurements  at  later 
times. 

We  remark  that  the  measured  values  used  in  Fig.  6  were  not  taken  from  Ref.  18.  While  our 
estimates  from  the  video  data  agreed  with  their  results  for  all  other  shot  cases,  our  estimates  for 
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Shot  11  were  25  to  30%  higher.  Since  the  digital  images  used  in  that  report  were  not  saved  and  hard 
copies  of  the  video  sequences  were  not  printed,  it  was  not  possible  to  reproduce  the  measurements 
presented  in  Ref.  18.  However,  we  attempt  to  justify  our  measurements  by  comparing  pictures  from 
the  video  sequence  and  computations  for  each  shot. 


0  0.2  0.4  0.6  0.8  1  1.2 

Time  (s) 


Fig.  5  —  Computed  and  measured  plume  heights  for  Shot  7 


Figure  7  shows  a  sequence  of  side-by-side  comparisons  of  video  frames  from  the  experiment  and 
computed  density  contours  from  Run  7-FI  of  the  plume  evolution  for  Shot  7.  At  early  times,  f  <  0.8, 
the  computations  under  predict  the  observed  plume  heights.  This  may  be  due  to  the  water  spalled 
upward  due  to  the  initial  shock  reflection.  The  outline  of  the  numerical  spray  region  appears  to 
provide  a  rough  approximation  to  the  outline  of  the  actual  plume.  The  emergence  of  the  secondary 
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plumes  through  the  spray  outline  coincides  at  t  =  0.8  in  both  the  computation  and  the  experiment. 
At  later  times,  the  secondary  plume  has  a  greater  horizontal  velocity  and  lower  vertical  velocity  in 
the  experiment  than  the  computation  indicates.  Still,  the  overall  qualitative  agreement  between  the 
computation  and  experiment  is  remarkable,  considering  the  complexity  and  duration  of  the 
phenomenon.  The  total  duration  of  the  plume  above  the  water  surface  is  slightly  under  predicted  by 
the  computations  as  indicated  at  t  =  5.0.  This  can  be  expected  due  to  the  drag  on  the  fine  water 
droplets  comprising  the  spray  in  the  plume  at  the  late  times. 

Figure  8  shows  a  sequence  of  side-by-side  comparisons  of  the  plume  evolution  for  Shot  11.  In  this 
case  the  early  time  plume  profiles  are  in  better  agreement  for  t  <  0.8  than  with  Shot  7.  Here,  the 
numerical  spray  region  more  closely  matches  the  opaque  spray  surrounding  the  pliune.  The  secondary 
plume,  which  can  be  seen  in  the  numerical  spray  region  at  t=  0.8,  has  not  yet  emerged  from  the 
spray  outline,  as  it  appears  in  the  video  frame.  Also,  as  with  Shot  7,  the  computed  secondary  plume 
rises  higher  than  in  the  experiment.  Since  this  secondary  plume  appears  to  be  comprised  primarily  of 
spray,  drag  may  be  effecting  its  motion  substantially.  The  overall  duration  of  the  plume  above  the 
surface  is  reproduced  well  by  the  computation,  as  indicated  at  r  =  4.4. 

Note  the  difference  in  scale  of  the  video  flumes  fl'om  the  two  shots  as  shown  in  Fig.  7  and  Fig.  8. 
The  grid  marks  in  each  case  were  based  on  matching  measured  plume  heights  at  corresponding  times 
from  video  fl’ames  taken  at  a  90°  angle  to  the  line  charge.  These  latter  frames  had  fiducials  so  that 
precise  scale  measurements  could  be  made  (see  Ref.  18).  As  an  additional  check  of  the  scale  used  in 
Figs.  7  and  8  observe  the  profile  of  the  trees  in  the  backgroimd.  The  field  of  view  (FOV)  for  Shot  1 1 
(Fig.  8)  is  obviously  narrower  than  for  Shot  7  (see  Fig.  7).  The  grid  scales  used  in  these  figures  closely 
match  the  difference  in  the  FOV.  When  this  difference  in  the  FOV  is  taken  into  aecount,  the  plume 
heights  for  Shot  11  agree  with  those  presented  in  Fig.  6  but,  as  previously  mentioned,  are  25  to  30% 
higher  than  reported  in  Ref.  18. 

Microwave  Data 

Measurements  of  the  plume  density  using  microwave  measurements  were  first  discussed  in 
Ref.  20.  These  measurements  were  based  on  the  amount  of  microwave  absorption  through  the  plume. 
Microwaves  were  sent  and  received  using  a  pair  of  3-ft  radius  parabolic  dishes  placed  on  either  side  of 
the  plume  at  equal  heights  above  the  water  surface.  To  compare  the  microwave  measmement  with 
the  computed  densities,  the  computed  values  were  integrated  within  a  cylindrical  region  of  radius  3  ft, 
corresponding  to  the  region  between  the  parabolic  dishes.  At  time  f,  this  integral,  representing  the 
total  amoimt  of  water  between  the  dishes,  can  be  expressed  as 

Iit%H,R,L)^r{H,R,L)=ridzf\dy\\p\x,y,z^^  (29) 

JH-K  •'-r(z)  J~L 

where  i?  =  3  is  the  radius  of  the  dishes,  H  is  the  height  of  the  center  of  the  dishes,  2L  is  the  distanee 
between  the  dishes,  and 


r(z)  = 


\4¥^^(z^liy  if\z-H\<R 
[  0  otherwise 


(30) 
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Fig.  8  —  Video  frames  and  computations  for  Shot  1 1 
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In  both  the  experiment  and  the  computations  the  value  for  L  was  sufficiently  large  to  contain  the 
entire  width  of  the  plume.  For  2-D  approximations,  p"  does  not  change  in  the  y'  direction  (parallel 
to  the  line  of  charges  and  perpendicular  to  the  line  between  the  microwave  dishes).  Therefore,  using 
symmetry  across  x  =  0,  it  follows  that 


H+R  L 

r{H,R,L)  =  4  J r(z)J p’'(x,y,z)dxdz. 

H~R  0 

This  integral  is  approximated  numerically  using  the  quadrature  formula 

nH,R,L)  -  II  =  ’ 

j=l 


(31) 


(32) 


where 


^  2 


{M)j  =  z]-z% 


(Ax),. 


zj  =  min(lT  +  R,  Zj^^ ),  z®  =  max(/f  -  R,  Zj ), 

=  max{  j:Zj<H-R}^,  +  ^}- 

This  mass  is  converted  into  an  equivalent  water  length  (EWL)  by 

nH,R,L) 


W{f,H,R,L)  =  W”{H,R,L)  =  - 


PqTcR^ 


(33) 


corresponding  to  the  length  of  water  filling  the  cylinder  having  an  equivalent  mass  as  the  plume 
intersected  with  the  cylinder.  The  cylinder  is  horizontal  with  its  axis  at  height  H,  with  radius  R,  and 
length  2L.  (For  the  microwave  data  the  length  of  the  cylinder  or  distance  between  the  dishes  is  not 
significant  since  this  distance  is  much  greater  than  the  width  of  the  plume.  That  is,  W  will  not  change 
if  L  is  sufficiently  large.) 

Figure  9  shows  the  computed  and  measured  microwave  data  for  Shot  7.  The  measwed  values  for 
this  shot  became  saturated  at  a  peak  value  of  1.56  as  indicated  by  the  flat  plateau  in  its  graph.  The 
graph  of  the  measured  data  begins  rising  approximately  0.15  s  earlier  than  the  computed  results  that 
corresponds  to  the  plume  height  data  presented  in  Fig.  5.  Since  the  plume  heights  fi-om  the  runs  with 
the  indented  surface  were  in  better  agreement  with  the  measurements,  only  data  from  those  runs  are 
presented  here.  At  early  times  (t<0.8),  only  the  central  plume  passes  through  the  microwave 
cylinder.  Both  computations  appear  to  severely  under  predict  the  amount  of  water  in  the  central 
plume  at  this  height.  Consider  the  structure  of  the  computed  plume  at  t=  0.5  in  Fig.  2.  The  density 
contours  indicate  a  thick  wall  of  water  rising  to  a  height  of  about  20  ft.  Above  this  thick  region  is  a 
much  thiimer  plume  resulting  from  the  initial  indentation.  According  to  the  computations,  the  thick 
part  of  the  plume  begins  to  thin  out  before  it  reaches  a  height  of  25  ft.  Therefore,  reducing  H 
generally  increases  W.  For  example.  Run  7-FI  yielded  JT(0.5,12.5,3.0,L)  —  3.3,  compared  to 
1F(0.5,25,3.0,Z,)  <  0.5,  as  indicated  in  Fig.  9.  The  peak  in  the  computed  values  at  approximately 
t  =  1.0  occur  as  the  secondary  plumes  pass  through  the  microwave  height.  This  peak  is  followed  by  a 
smaller  peak  as  the  water  at  the  top  of  the  secondary  plume  falls  back  downward.  Since  the  secondary 
plume  was  not  ejected  upward  as  high  in  the  coarse  grid  run,  the  smaller  peak  appears  substantially 
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earlier  (^  =  1.3)  for  Run  7-CI  than  the  time  (t  =  2.1)  it  appears  for  Run  7-FI.  At  later  times  {t>  3), 
the  fine  grid  run  and  the  measured  data  are  in  better  agreement. 


0  1  2  3  4  5  6 

Time  (s) 

Fig.  9  —  Computed  and  measured  microwave  data  for  Shot  7 


Figure  10  displays  the  microwave  results  for  Shot  11.  Here  the  initial  rise  in  the  graphs  for  f  <  0.5 
are  all  in  agreement.  However,  while  the  measured  data  rises  and  then  falls  (consistent  with  the  peak 
computed  plume  heights  shown  in  Fig.  4),  the  computed  data  for  Run  11 -FI  peaks  at  1.0  as  the 
secondaiy  plume  rises  above  25  ft  (cf..  Fig.  8).  While  the  computed  equivalent  water  length  values 
range  between  1.1  and  1.7  for  1  <f  <2.1,  the  measured  data  range  between  2.0  and  3.2  during  the 
same  time  period. 


Fig.  10  —  Computed  and  measured  microwave  data  for  Shot  11 


Probe  Data 

In  addition  to  the  microwave  data,  plume  densities  were  measured  using  conductivity  probes. 
These  probes  were  developed  by  Phillips  and  Scott  [21]  and  consist  of  two  parallel  stainless  steel  rods 
whose  conductivity  is  linearly  proportional  to  its  imwetted  length.  For  these  tests,  fifteen  probes 
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were  suspended  on  a  cable  perpendicular  to  the  line  of  charges  at  a  nominal  height  of  12.5  ft.  The  use 
this  technology  for  plume  density  measurements  was  conceived  of  and  previously  used  by  Lipton  in 
Ref.  22.  However,  the  measured  height  above  the  charges  was  only  11.5  ft  as  described  in  Ref.  23. 
The  probes  were  spaced  1  foot  apart  with  a  central  probe  directly  above  the  line  of  charges. 

Making  comparisons  of  density  measurements  at  the  individual  locations  is  not  meaningful  due  to 
the  instabilities  inherent  in  the  plume  formation.  Indeed,  there  is  very  little  agreement  between  probe 
pairs  equally  distant  on  either  side  of  the  charges.  However,  some  success  was  achieved  by  making 
comparisons  using  an  integral  norm.  As  with  the  microwave  comparisons,  we  define  the  effective 
water  length  W  using  (29)  and  (33),  taking  the  limit  as  i?  0.  This  yields 

1  ^ 

W(t,HAL)  =  —  J  P(f,x,yp,H)dx,  (34) 

Po  -L 


where  H  is  the  height  of  the  probe  line,  jp  the  horizontal  offset  from  the  central  charge  (or  center  of 
the  line  charge),  and  2L  is  the  distance  between  the  first  and  last  probe  in  the  line.  Since  the  probes 
are  located  at  discrete  points,  this  quantity  is  approximated  using  the  trapezoid  rule  quadrature 


N~\ 

Pi  /2+^p,. +  PAf  /2  , 


i=2 


(35) 


where  S  is  the  uniform  spacing  between  the  N  probes  and  p.  is  the  density  at  the  ith  probe  location. 
The  use  of  (36)  with  5=  1  and  A  =  15,  corresponding  to  the  actual  probe  locations,  is  referred  to  as 
“point  line  integration”  (PLI).  Later,  we  will  also  consider  an  approximation  to  Eq.  (34)  based  on 
values  for  the  densities  inside  each  computational  cell,  with  values  for  L  greater  than  the  plume 
widths.  This  will  be  referred  to  as  “full  line  integration”  (FLI). 

The  computational  data  was  integrated  at  i/=  12.5,  as  opposed  to  the  11.5  ft  height  of  the  cable. 
However,  the  cable  was  moved  by  the  plume  during  the  experiment  so  precisely  that  inatching  the 
height  of  the  probe  cable  was  not  feasible.  Figure  1 1  shows  a  comparison  of  the  equivalent  water 
lengths  Ws,  using  Eq.  (35)  on  the  probe  data  and  the  computational  data.  While  the  peak  values  of  Ws 
agree  to  within  15%,  the  peak  of  the  computed  values  occurs  at  approximately  t  =  0.5  s,  while  the 
measured  peak  occurs  at  approximately  t  =  1.2  s.  Compared  to  the  microwave  data,  both  the 
computed  and  measured  probe  data  indicate  much  lower  fV  values  for  t>  2.  This  is  because  after  this 
time  most  of  the  plume  has  spread  out  beyond  the  15  ft  width  straddled  by  the  probes.  Furthermore, 
the  probes  open  downward  and  are  not  able  to  detect  water  falling  downward  due  to  a  “shadowing” 

effect. 

Figure  12  displays  the  comparison  for  Shot  1 1 .  Here,  the  computations  are  substantially  higher 
than  the  measurements  at  almost  all  times.  As  the  water  in  the  plume  is  broken  up  into  droplets  and 
spray,  the  behavior  of  the  probes  is  not  well-knovra.  In  general  droplets  smaller  than  the  width  of  the 
fork  in  the  probes  will  not  be  detected.  This  may  partially  explain  the  large  discrepancy  in  the  data. 

The  integrity  of  all  the  data  can  be  checked  by  comparing  results  from  each  measurement  for 
one  of  the  shots.  Fig.  13  shows  a  comparison  of  the  data  for  Shot  7  for  t<l.  For  t  <  0.8,  the  plume 
is  comprised  of  only  the  central  water  column,  surrounded  by  spray  from  the  initial  shock  reflection 
and  Rayleigh-Taylor  instability  at  the  free  surface.  Therefore,  it  can  be  expected  that  the  density 
probes  span  most  of  the  water  in  the  plume  at  these  times.  However,  as  seen  in  Fig.  13,  the 
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integrated  measured  probe  data  is  not  only  far  below  the  computed  data  (Run  7-FI  (PLI)  and  (FLI)), 
but  it  is  also  well  below  the  measured  microwave  data  recorded  at  over  twice  the  height  of  the  probe 
line.  That  the  probes  are  spaced  sufficiently  close  to  resolve  the  plume  structure  and  their  14-ft  span 
is  sufficiently  wide  during  these  times  is  supported  by  the  similarity  of  the  computational  data  using 
either  point  line  integration  at  only  the  probe  locations  (PLI)  or  full  line  (FLI)  integration  at  all  cell 
locations.  These  two  integration  formulas  begin  to  diverge  after  t  =  0.8  as  the  wider  secondary  plume 
reaches  the  probe  height.  Another  disturbing  feature  in  Fig.  13  is  that  the  microwave  data,  measured 
at  25  ft  begins  to  rise  at  about  the  same  time  as  the  probe  data,  measured  at  11.5  ft.  These 
discrepancies  suggest  that  inaccuracies  in  the  measurements  may  be  as  significant  as  with  the 
computations. 


Fig.  12  —  Computed  and  measured  probe  data  for  Shot  11 
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Optimal  Depth  Study 

Despite  the  relatively  poor  agreement  between  the  computations,  microwave  data,  and  probe 
data,  an  optimality  study,  based  on  the  computations,  is  presented.  The  validity  of  this  study  is  based 
on  the  agreement  in  the  measured  plume  heights  and  the  qualitative  agreement  between  the  video 
frames  and  the  computed  density  contours.  The  use  of  microwave  and  probe  measurements  has  not 
been  validated  for  plume  density  measurements  independently.  Therefore,  even  though  the 
computations  may  (or  may  not)  disagree  with  the  actual  plume  water  content  at  any  individual 
charge  depth,  the  dependence  of  the  plume  density  on  the  charge  depth  may  still  be  accurately 
predicted. 


The  computations  were  performed  with  initial  charge  depths  between  1  and  21  ft.  In  particular 
runs  were  made  at  depths,  1,  1.5,.. .,4.5,  5.0,  5.2,.. .,9.8,  10.0,  11.0,...,15.0,  17.0,  19.0,  and  21.0. 
The  initial  conditions  used  were  the  same  as  those  based  on  the  empirical  laws  described  previously, 
modeling  a  line  charge  of  M=1.25  Ibs/ft  of  C-4.  Distances  (including  charge  depths)  may  be 
nondimensionalized  using  the  free-field  maximum  bubble  radius.  However,  the  effects  of  gravity  do 
not  scale,  so  the  “nondimensionalized”  results  presented  here  can  only  be  expected  to  be  valid  for 
charges  producing  bubble  energies  within  an  order  of  magnitude  of  these  C-4  charges.  Initially,  the 
computations  used  an  indented  free  surface  with  the  radius  of  the  indentation  given  by  if/=1.03J, 

whenever  the  scaled  depth  is  less  than  one,  that  is,  C  =  -^<1.  As  with  the  computations  for  Shot 


/111 


7  and  Shot  1 1 ,  the  runs  were  started  on  a  grid  that  had  a  uniform  fine  grid  region  of  cells  Avith  size 
^  =  0.05,  surrounding  the  initial  charge  location.  In  all  cases,  the  uniform  region  in  the  horizontal 
direction  was  restricted  to  the  vertical  strip,  0<x<2.  The  uniform  region  in  the  vertical  direction 
depended  on  the  initial  charge  depth,  as  indicated  in  Table  3.  In  the  cases  when  the  top  of  the 
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uniform  region  was  below  z  =  l,  cells  of  height  2/i,  were  added  to  extend  the  computational  domain 
to  z  =  l.  As  before,  these  initial  grids  were  used  for  0< ?</■,.  =0.007,  after  which  the  solution  was 
remapped  onto  the  second  grid.  The  uniform  grid  region  for  the  second  grid  had  cell  size  =  0.4, 
corresponding  to  the  “coarse”  grid  used  previously.  The  imiform  grid  region  used  for  the  second  grid 
was  the  vertical  strip  0  <  x  <  12  intersected  with  the  horizontal  strip  indicated  in  Table  3. 


Table  3  —  Uniform  Grid  Regions  for  Depth  Dependence  Runs 


Initial  Grid 

Second  Grid 

d<5 

5<d<l0 

n<d<l5 

d>l5 

-7<z<l 
-12<z<-3 
-17<z<-9 
-d-2 <  z  <-d+2 

-21<z<29 

-21<z<29 

-25<z<29 

-d-\2<z<29 

In  this  study  the  expression  “optimal”  will  refer  to  some  functional  F  of  the  equivalent  water 
length  fV  as  a  function  of  the  depth,  d.  Ideally,  we  seek  dop„  which  satisfies 

F{W{t,H,R,L,d,^))=  m^F{}V{t,Ii,R,L,d)).  (36) 


Two  specific  forms  for  F  will  be  used.  The  first  is  simply  the  time  integral  of  W,  namely, 

which  has  units  of  mass-time.  The  second  functional  measures  the  length  of  time  that  W  is  greater 
than  or  equal  to  some  threshold.  More  precisely, 

Ud,E^J)  =  F^mt,Hj,,R,,,oo,d\E^,T)  =  m{{t-.W{t,H,,,R,,,oo,d)>E^}n{0<t<T}),  (38) 

where  m  refers  to  the  Lebesgue  measure  (length)  of  the  set.  The  threshold  value  E^y  may  relate  to 
some  predetermined  value  for  which  the  plume  makes  an  effective  barrier. 

Three  choices  for  the  pair  (H,  R)  are  (12.5,  0.5),  (16.5,  5.5),  and  (25.0,  0.5).  These  will  be 
referred  to  as  low,  average,  and  high,  respectively.  The  “low”  and  “high”  choices  correspond  to  what 
would  be  encoimtered  by  a  one-foot  diameter  missile  at  heights  of  12.5  and  25  ft,  respectively.  The 
“average”  represents  a  weighted  average  between  heights  11  and  22  ft.  This  roughly  corresponds  to 
the  scaled  heights  between  1  and  2  A„ax- 

A  graph  of  the  function  ftid,T)  for  the  three  choices  of  (H,  R),  is  shown  in  Fig.  14  below,  with 
T  =  6.0  (the  entire  plume  duration).  Due  to  the  oscillations  in  the  computed  data,  an  approximation 
of  dop,  cannot  be  confidently  determined.  However,  some  of  the  features  of  these  graphs  become 
clearer  after  tbe  data  has  been  smoothed.  The  smoothing  of  the  data  was  performed  using 
g^  =  ^'“(g),  where  A{g\  =  (g,._,  +  4g,  +  g,+i)/6  and  g  is  the  vector  of  values  of fi{d,T)  interpolated  (if 
necessary)  to  the  uniform  grid  t/,  =1  +  0.2/  for  /  =  0,...,100.  The  values  at  the  endpoints  /  =  0  and 
/=  100  are  held  fixed  by  the  smoothing. 
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Fig.  14  —  Smoothed  and  raw  computed  values  forfi{d,T) 


A  phenomenological  argument  explaining  all  of  the  relative  extrema  shown  in  the  data  cannot  be 
made  at  this  time.  However,  the  three  clear  relative  maxima,  appearing  in  the  “average”  case,  can  be 
partially  explained  by  a  careful  examination  of  the  computed  results.  At  shallow  depths,  C  <  0.5 
{d  <  6),  the  bubble  can  be  expected  to  vent  before  its  first  period  (time  of  minimum  volume).  For 
very  shallow  depths  {d  <  2.0),  the  venting  occurs  early  and  the  amount  of  water  ejected  upward  is 
almost  exactly  proportional  to  amount  initially  above  the  bubble  (charge).  As  the  depth  increases, 
the  time  of  venting  also  increases.  Eventually,  venting  will  not  occur  until  after  the  bubble  attains  a 
pressure  below  the  ambient  pressure.  During  this  time,  the  water  above  the  bubble  is  accelerated 
downward,  thereby  decreasing  the  amount  ejected  upwards.  This  explains  the  first  relative  maximum 
at  if  =  2.5  (C  =  0.21).  When  venting  is  delayed  until  after  the  first  bubble  maximum  volume,  the 
velocities  of  the  surrounding  water  point  toward  the  bubble’s  center  causing  it  to  contract.  A  high 
pressure  region  forms  above  the  bubble  due  to  the  water  rushing  inward  in  the  relatively  thin  layer 
between  the  bubble  and  the  air.  This  causes  the  formation  of  the  central  plume  ((cf..  Figs.  7  and  8) 
t  =  0.4  and  0.5).  Due  to  the  relatively  thin  layer  of  water  above  the  bubble  at  this  time  and  Rayleigh- 
Taylor  instability,  some  “fingering”  of  the  free  surface  can  pierce  the  bubble  causing  some  venting  to 
occur.  The  amount  of  water  in  the  plume  depends  on  how  late  venting  occurs.  At  depths  d  >  6.0 
(C  >  0.5),  the  bubble  no  longer  vents  during  its  first  pulse.  As  the  depth  increases,  the  central  plume 
thickens  but  attains  lower  maximum  heights.  The  second  relative  maxima  in  /;  for  the  “average 
choice  occurs  at  approximately  d  =  7.6,  (C  =  0.66).  Adding  to  the  total  water  ejected  are  the 
secondary  plumes  that  appear  during  the  bubble’s  second  expansion  (cf.,  Figs.  7  and  8,  t  =  0.8).  While 
the  total  amount  of  water  in  the  central  plume  appears  to  diminish  after  d>  8.0,  the  amounts  ejected 
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by  the  secondary  and  tertiary  pulses  increase,  and  the  global  maximum  appears  at  d=  14.0,  (C  =  1.3) 
for  the  “average”  choice.  At  the  “low”  height,  the  global  maximum  occurs  at  d=l5  and,  at  the 
“high”  height,  it  occurs  at  the  shallower  depth  of  13.  This  can  be  explained  by  tertiary  plumes 
contributing  more  water  and  rising  to  maximum  heights  between  12.5  and  22  ft. 

The  theory  described  above  is  also  supported  by  the  graph  of  the  smoothed  second  functional  f2 
defined  by  Eq.  (38).  Figures  15  and  16  display  graphs  of  /2(J,1.5,6.0)and  /jCrf,!. 5,2.0),  respectively. 
The  similarity  in  these  two  sets  of  graphs  for  d<  10  demonstrates  that  for  these  depths  most  of  the 
water  appears  in  the  plumes  for  the  first  2  s  after  the  detonation.  However,  for  12,  a  substantial 
amount  of  water  is  ejected  upward  after  2  s.  This  corresponds  to  times  after  several  bubble  pulses 
have  occurred.  The  relative  maximum  at  (i=18.8  at  the  “low”  height  for  1.5, 6.0)  is  almost 

twice  its  corresponding  value  for  ^(i/,1. 5,2.0).  This  indicates  that,  for  depths  near  d=l9,  a  plume 
rising  just  above  12.5  ft  is  ejected  upward  shortly  before  t  =  2  and  contains  a  substantial  amount  of 
water. 
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Fig.  15  —  Smoothed  computed  values  for f2{d, l.5,6.0) 
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Fig.  16  —  Smoothed  computed  values  for f2{d,l. 5,2.0) 
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Validations  of  the  Three-Dimensional  Model 

The  BUB3D  code  was  used  to  model  the  effects  of  using  discrete  charges  as  opposed  to  the 
continuous  line  charge  approximated  with  the  2-D  model.  The  goal  of  this  computational 
experiment  was  to  predict  the  difference  in  the  plume  structure  observed  between  the  use  of  discrete 
or  continuous  line  charges.  This  difference  is  conjectured  to  be  caused  from  shock  interactions  at  the 
free  surface  between  the  discrete  charges  shortly  after  a  (nearly)  simultaneous  detonation.  The 
incompressible  liquid  model  was  initialized  using  the  empirical  shock  model  depicted  in  Fig.  1.  In 

particular,  the  choices  Ri  =  1.03^/  and  Rp  =  -^{S Ilf  were  used  to  initialize  the  density. 

The  assumption  that  the  line  charge  is  of  infinite  extent  was  retained  from  the  2-D  model,  so 
that  only  one  variable  (initial  charge  shape)  was  changed.  This  approximation  was  implemented, 
using  an  extra  symmetry  plane,  centered  between  two  discrete  charges.  A  second  symmetry  plane, 
parallel  to  the  first  and  cutting  through  the  center  of  the  initial  charge  location,  was  also  used. 
Finally,  the  model  includes  a  third  symmetry  plane,  containing  the  line  of  charges.  This  third  plane 
corresponds  to  the  symmetry  plane  used  in  the  two  dimensional  computations.  If  an  initial  charge  is 
located  at  coordinates  (x,  y,  z)  =  (0,  0,  -d)  and  the  initial  standoff  distance  between  the  charges  is  S', 
then  the  three  symmetry  planes  described  above  are  located  aty  =  5/2,  y  =  0,  and  x  —  0,  respectively. 

Using  the  conditions  for  Shots  2,  6,  and  9,  namely,  d—  8.2,  10  lbs  of  Composition  C-4,  it 

follows  from  Eqs.  (18-23)  that  4,, „  =0.61617,  4n^  =  9.3,  and  =1535.7P^.  (Note,  this  is  the  same 
initial  pressure  as  for  the  2-D  model  for  Shot  7.)  As  with  the  2-D  model,  a  two-grid  solution 
procedure  was  used.  The  initial  grid  contained  a  fine  uniform  region  of  cell  cubes  of  size  —  0.1  in 
the  region  0<x<2,  intersected  with  0<y  <4  and  -10.2 < z <~6.2.  Above  this  region  the  cells  were 
stretched  vertically  to  a  maximum  size  of  =0.2,  using  46  cells,  extending  the  domain  to  z  =  2.0. 
The  grid  was  stretched  downward  to  z  =  —50,  using  an  additional  12  cells.  In  the  x-direction,  cells  were 
stretched  to  x  =  50,  using  an  additional  10  cells.  Overall,  the  initial  grid  was  comprised  of 
30  X  40  X  98  cells. 

The  second  grid  had  a  resolution  corresponding  to  the  “coarse”  grids  used  in  the  2-D  study.  The 
uniform  region  was  composed  of  cubic  cells  of  size  —  in  the  region  0  <  x  <  12  intersected  with 
0  <  y  <  4  and  -20  <  z  <  30 .  Additional  stretched  cells  extended  the  domain  down  to  z  =  -100,  up  to 
z  =  150,  and  across  to  x  =  100.  The  second  grid  used  a  total  of  50  x  10  x  185  cells. 

Figure  17  shows  a  composite  of  video  images  from  Shot  6  and  the  3-D  computations.  The  view 
displayed  in  these  images  is  the  “front”  view,  perpendicular  to  the  line  of  charges  and  the  images 
shown  in  Figs.  7  and  8.  Density  plots  from  the  computations  overlay  the  video  images  in  the  range 
0<y  <28.  These  plots  were  created  by  reflecting  and  copying  plots  from  the  actual  computational 

region  (0<y<4)by  symmetry. 

The  video  image  at  time  t  =  0  corresponds  to  a  time  shortly  after  the  charge  detonations.  The 
density  plot  overlaying  this  image  corresponds  to  the  computation  at  ^  =  0.0035,  the  initial  time  for 
the  remapped  solution  on  the  second  grid.  At  this  time  the  bubbles  are  still  discrete  and  nearly 
spherical.  The  indentations  at  z  =  0  for  the  empirical  model  of  the  shock  interactions  can  be 
observed.  The  shock  interactions  may  be  clearly  observed  as  white  “bumps”  above  the  line  of  charges 
at  the  midpoints  between  the  charges.  Curved  white  lines  emanating  from  these  points  on  either  side 
of  the  charge  line  are  also  observable  in  this  image.  The  curve  in  these  lines  is  due  to  slightly 
asynchronous  detonation  (from  right  to  left)  of  the  charges.  At  ^  =  0.1,  seven  distinct  finger 
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plumes”  can  be  seen  rising  to  a  height  of  about  20  ft,  compared  to  15  ft  in  the  computations.  The 
computed  shorter  plumes  directly  above  the  initial  charge  locations  arise  from  the  initial  indentation 
of  radius  In  the  video  it  is  difficult  to  see  this  detail  in  the  plume  structure.  The  computed  bubble  is 
now  cylindrical  shaped  separated  by  a  thin  membrane  at  the  symmetry  wall.  The  bubble  expands 
attaining  a  maximum  volume  at  t  =  0.266,  and  the  contraction  can  be  clearly  observed  at  t  =  0.4.  The 
bubble  continues  its  collapse  imtil  shortly  after  t  =  0.5.  The  height  of  the  plume  fingers  is  under 
predicted  by  the  computations  for  t<0.3  and  slightly  over  predicted  for  t>0.3.  At  the  later  times, 
the  plume  velocity  can  be  expected  to  be  influenced  by  air  resistance.  Air  resistance  is  not  modeled  in 
the  computations.  This  is  a  likely  cause  of  the  over  production  of  the  plume  heights  for  t  >  0.3. 


1=0.3  1=0.4  t=0.5 


Fig.  17  —  Composite  of  video  images  and  computed  density  plots  for  Shot  6 


A  comparison  of  the  measured  and  computed  plume  heights  is  shown  in  Fig.  18.  Also  included 
were  measured  plume  heights  from  Shots  2  and  9  (same  depth  and  standoff  as  Shot  6  but  using  only  5 
charges).  The  plume  heights  from  Shots  2  and  9  are  slightly  lower  than  for  Shot  6,  and  they  are  in 
better  agreement  with  the  computation  for  t  <  0.3  and  in  worse  agreement  for  t  >  0.3. 

Figure  19  displays  perspective  images  of  the  3-D  computation.  In  these  frames,  two  isocontours 
of  density  are  displayed.  For  the  plume  above  z  =  0,  the  isosurface  p  =  0.001pQ  is  rendered  in  light 
gray,  while  for  the  bubble  below  z  =  0  the  isosurface  p  =  O.Sp^  is  rendered  in  dark  gray.  At  t  =  0.2,  the 
bubble  is  clearly  seen  to  be  cylindrical  in  shape.  A  2-D  jet  can  be  seen  piercing  the  top  of  the 
cylindrical  bubble  at  the  start  of  the  collapse  phase  at  t  =  0.4.  The  ejection  of  the  secondary  plumes 
can  be  seen  at  times  t  =  1.0  and  1.2.  These  secondary  plumes  appear  later  in  this  calculation  than 
in  the  analogous  2-D  computation  (cf.,  Fig.  7)  due  to  the  bubble  venting  into  the  air  region  during  its 
collapse.  When  the  bubble  reformed  it  continued  to  collapse  but  had  less  energy  for  the  second 
expansion. 
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Fig.  19  —  Computed  bubble  and  plume  isocontours  for  Shot  6 
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SUMMARY  AND  CONCLUSIONS 

This  report  presented  comparisons  between  computational  and  experimental  measurements  of 
plumes  produced  by  underwater  explosions.  Of  particular  interest  was  the  study  of  charge 
configurations,  which  create  a  plume  “barrier”  by  ejecting  a  “wall  of  water”  above  the  surface.  Such 
an  effect  occurs  after  either  the  nearly  simultaneous  detonation  of  discrete  charges  placed  sufficiently 
close  together  in  a  line,  or  the  detonation  of  a  continuous  line  of  charges. 

The  overall  hydrodynamic  phenomena  involved  in  this  process  is  extremely  complicated.  The 
bubble  formed  from  underwater  detonation  can  undergo  several  cycles  or  periods  before  it  vents. 
During  each  cycle,  hydrodynamic  instabilities  can  be  expected  at  both  the  bubble-water  and  air-water 
interfaces.  Rayleigh-Taylor  instabilities  occur  whenever  a  denser  material  (water)  is  accelerated  into  a 
less  dense  material  (bubble  gas,  or  air).  These  instabilities  occur  at  the  bubble  interface  when  the 
bubble  is  near  its  minimum  volume  and  the  bubble  pressure  is  above  the  ambient  hydrostatic  pressure. 
They  also  occur  at  the  air-water  interface  as  the  bubble  begins  its  contraction.  Such  instabilities  make 
pointwise  density  comparisons  meaningless.  Combined  with  the  long-time  behavior  of  the  overall 
dynamics,  the  numerical  accuracy  of  the  computations  is  also  an  issue.  Nevertheless,  the  overall 
dynamics  of  the  plume  formation  and  secondary  plume  ejection  is  reproducible  experimentally,  and 
as  was  demonstrated  in  this  paper,  computationally  as  well. 

Numerical  accuracy  of  the  computations  was  studied  by  comparing  the  predicted  d)mamics  using 
two  different  grids  (fine  and  coarse),  with  sizes  differing  by  a  factor  of  two.  For  relatively  early  times 
{t  <  1 .0)  the  differences  in  the  solutions  were  very  small  (cf..  Fig.  2)  hut  became  more  apparent  at 
later  times.  These  differences  were  also  studied  quantitatively  by  comparing  computed  plume  heights 
(cf..  Figs.  3  and  4).  Comparisons  were  made  with  and  without  an  indented  free  surface  used  to 
empirically  model  the  shock  effects. 

The  use  of  an  empirical  model  (cf..  Fig.  1)  for  shock  effects  at  the  air-water  surface  was  found  to 
be  critical  for  accurate  predictions  of  the  plume  heights.  For  continuous  line  charges,  the  model 
accounts  for  the  reflection  of  the  cylindrical  shock  off  the  air-water  interface  by  initializing  the 
interface  with  a  small  channel  in  the  water  directly  above  the  line  charge.  When  this  model  was  used, 
the  predicted  pliune  heights  were  initially  lower  than  those  observed  (for  r<0.3)  probably  due  to 
water  spalled  upward  from  the  shock  reflection  unaccounted  for  in  the  model.  At  later  times  the 
computed  plume  heights  eventually  exceed  the  observed  heights  (cf..  Figs.  5  and  6).  This  over 
prediction  is  conjectured  to  be  caused  by  the  lack  of  air  resistance  in  the  computational  model. 

By  also  including  the  effects  of  shock  interaction  between  discrete  charges,  the  plume  heights  for 
a  discrete  line  charge  case  accurately  matched  the  observed  heights  using  a  3-D  computational  model 
(see  Fig.  18).  In  the  discrete  charge  case,  the  plume  heights  are  generally  higher  than  with  an 
equivalent  line  charge  due  to  the  “fingering”  of  the  plumes  between  the  charge  locations.  This 
fingering  effect  was  reproduced  by  the  3-D  computation  using  the  empirical  shock  model  (see 
Fig.  17). 

Secondary  plumes  are  ejected  (for  scaled  charge  depths  C>0.5)  for  line  charges  and  single 
discrete  charges.  For  discrete  charges,  secondary  (radial)  plumes  were  first  accurately  predicted  in 
Ref.  2.  For  continuous  line  charges,  the  emergence  of  secondary  plumes  has  been  observed  and 
accurately  predicted  in  this  study  (see  Figs.  7  and  8).  These  secondary  plumes  were  also  predicted 
using  the  3-D  model  (see  Fig.  19).  This  computation  also  demonstrated  the  cylindrical  shape  of  the 
merged  bubbles  during  their  first  expansion. 
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Comparisons  of  computed  and  measured  plume  densities  were  unfortunately  less  successful.  The 
computations  predicted  significantly  less  density  in  the  plume  than  the  amounts  determined  using  the 
microwave  absorption  measurements  (see  Figs.  9  and  10).  However,  the  computations  predicted 
either  more  (see  Fig.  12)  or  roughly  the  same  (see  Fig.  11)  total  amoxmt  of  water  in  the  plume  as 
determined  by  the  conductivity  probe  measurements.  The  disagreement  of  the  measured  results,  due 
to  the  measurement  techniques  and  the  lack  of  calibration,  undermines  their  use  as  a  validation  tool 
for  the  computations. 

The  optimal  depth  study  was  complicated  by  the  relative  unsmoothness  of  the  results  as  a 
function  of  depth  (cf..  Fig.  14).  Even  when  the  results  were  regularized,  several  relative  maxima 
appeared.  Furthermore,  the  results  were  dependent  on  both  the  measure  used  and  the  height  at  which 
the  densities  were  integrated.  For  the  range  of  depths  studied  (0.1  <  C<2.2),  three  relative  maxima 
appeared  consistently  for  all  the  cases  considered.  Conjectures  were  made  to  explain  these  relative 
maxima  but  were  not  verified  by  either  additional  computational  or  experimental  tests. 

Improvements  in  the  predictive  model  can  be  made  in  several  ways.  One  improvement  would  be 
to  eliminate  the  empirical  conditions  used  for  the  initial  conditions,  using  a  validated  code  which 
includes  compressible  effects  for  the  water  region.  The  compressible  code  would  need  to  be  run  for 
only  a  relatively  short  duration  until  the  shock  and  reflected  rarefaction  wave  are  sufficiently  far 
from  the  initial  charge  location.  After  this  time,  the  results  of  the  compressible  code  could  be  used  as 
initial  data  for  the  currently  used  BUB2D  or  BUB3D  codes.  Another  relatively  straightfonvard 
improvement  would  be  the  inclusion  of  a  model  for  the  “air”.  In  the  current  model,  the  air  is  simply 
a  void  with  uniform  density.  Treating  the  air  as  a  second  incompressible  species,  with  a  prescribed 
density,  may  yield  more  accurate  late-time  plume  feature  predictions. 
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